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ABSTRACT 

We present a new flexible estimator for comparing theoretical templates for the pre- 
dicted bispectrum of the CMB anisotropy to observations. This estimator, based on 
binning in harmonic space, generalizes the "optimal" estimator of Komatsu, Spergel, 
and Wandelt by allowing an adjustable weighting scheme for masking possible fore- 
ground and other contaminants and detecting particular noteworthy features in the 
bispectrum. The utility of this estimator is illustrated by demonstrating how acous- 
tic oscillations in the bispectrum and other details of the bispectral shape could be 
detected in the future PLANCK data provided that Jnl is sufficiently large. The char- 
acter and statistical weight of the acoustic oscillations and the decay tail are described 
in detail. 
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1 INTRODUCTION 

At present tlie primordial blackbody component of the 
cosmic microwave background (CMB) appears very nearly 
Gaussian, as one would expect from inflationary models for 
the very early universe. To a first approximation, because of 
the weak coupling of the physics involved, inflation predicts 
a nearly scale-invariant spectrum of primordial density per- 
turbations whose imprint on the CMB is completely charac- 
terized by the power spectrum Cf^, where I is the multipole 
number and (A, B = T, E), according to the Gaussian dis- 
tribution: 



n(2vr)-'/'det-^/' [Ci] 



X exp 



1 T ^-1/2 



where 



CLl,m — \ E 



C, 



TT (jTE 



(1) 



(2) 



with T denoting temperature and E the i5-mode polariza- 
tion. 

It has however been noted that even in single scalar field 
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inflationary models small nonlinear corrections do arise, at 
lowest order in the bispectrum (Maldacena' '2003 Acqua- 



viva et al. 20031. In models with several scalar fields even 



larger and potentially detectable degrees of non-Gaussianity 
are possible ([Bartolo, M atarrese fc Riotto 20021 'Rigopot^j 
los, Shellard fc Van Tent 2 006 2007; Naruko fc Sa saki 2009) 
Byrnes et al.|2008[ |2009|fHuang|2009| ! 



Non-Gaussianity manifests itself in odd n-point correla- 
tion functions or in the connected even n-point correlation 
functions, from which the trivial part expressible as com- 
binations of two-point correlation functions has been sub- 
tracted away. The extent of departures from Gaussianity 
can be characterized by ratios of higher-order correlation 
functions and the appropriate combination of two-point cor- 
relation functions ( Bernardeau et aLl |2002). The evolution, 
for the most part linear, of the primordial fluctuations of 
the inflaton field, involving both gravity and hydrodynam- 
ics, leads to CMB anisotropies whose statistical properties 
mirror those of the primordial fluctuations. Consequently, 
by studying higher-order correlation functions of the CMB 
anisotropies, we can detect and characterize any primordial 
non- G aussianity. 

The lowest order such statistic is the bispectrum, or 
three-point correlation function in Fourier space. The bis- 
pectrum has been shown to be an optimal statistic for 
measuring non-Gaussianity in the sense that the signal-to- 
noise squared of the non-Gaussianity estimator based on the 
three-point correlation function dominates over all higher- 
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order estimators ( Babich||2005 1. Consequently it would also 
be significantly easier to constrain. Assuming that the char- 
acteristic amplitude of the bispectrum is much larger than 
that of the trispectrum, or the four-point correlation func- 
tion in Fourier space, the non-Gaussianity of the density 
fluctuations from inflation for the single-fleld case may be 
approximated by the following "local" ansatz for the grav- 



Matarrese fc Riotto|[2004l ICreminelU fc Zaldarri"aga||2004a| 



Nitta et al. 



20091, as well as higher orders in the gravi- 



tational potential during recombination (Bartolo & Riotto 
|2008( ), can generate non-Gaussianity. Other sources of non- 
primordial non-Gaussianity include secondary anisotropies 
such as weak lensing via the cross-correlation with the un- 
lensed CMB arising from the integrated Sachs- Wolf effect 



20011 



itational potential (Verde et al. 2000 Komatsu & Spergel (Smith & Zaldarriaga 2006; Hanson et al.|[2009| Mangilli & 



$(x) = <1>g(x) + /ml [$g'(x) - {$g'(x))] 



(3) 



where the nonlinearity parameter /nl may be considered 
small. While the three-point correlator of the linear or Gaus- 
sian part <&G is zero (because in the quantum picture it in- 
volves the expectation value of a product of three creation 
and annihilation operators), this is no longer true when the 
quadratic correction is taken into account. The bispectrum 
will be proportional to /nl, which determines the size of 
non-Gaussianity. While this form is only approximate for 
the single-field case, for the observationally more promising 
multi-field models for which potentially observable values of 
/nl can be generated, this ansatz is very accurate because 
the bulk of the non-Gaussianity is imprinted after horizon 
crossing. 

The primordial bispectrum is characterized by an over- 
all amplitude /ml and the shape. Translational and rota- 
tional invariance reduce the bispectrum to a function of the 
lengths only of the three wave vectors. A classification of 
bispectra shapes has been proposed by |Fergusson fc Shel-| 
|lard| [2009) . Local and warm models peak on squeezed trian- 
gles, the equilateral models on equilateral triangles, and the 
folded/fiat models on fiattened triangles. Models in which 
perturbations are generated outside the horizon produce a 
bispectrum that peaks on squeezed triangle configurations, 



examples of which are multi-field infiation ( [Lyth, Ungarelli 
fc Wands|2003l [Bartolo, Matarrese fc Riotto|2004||Rigopo'i> 
los, SheUard fc Van Tent||2007|) . [For" a review, see IBartolo, 



Komatsu, Matarrese fc Riotto| ( |2004[ ).] The bispectrum of 



standard single-field inflation can be viewed as a superpo- 
sition of the local shape and the equilateral shape, both 



terms, however, being slow-roll suppressed ( 


Maldacena 


2003 


Fergusson & Shellard||2009| |Senatore, Smith fc Zaldarriaga] 


2009 


I. In particular, for the squeezed triangle configurations 



in single-field inflation the non-Gaussian signal would be 
proportional to the tilt of the power spectrum and thus im- 
ply a strong deviation from scale invariance ( [Creminelli fc| 
Zaldarriaga|2004b Cheung et al.|2008 1 . A bispectrum peak- 



ing on equilateral triangles is predicted in models with non- 
standard kinetic terms in the infiaton Lagrangian, which is 
the case of multi-field inflation models such as DBI inflation 
(Alishahiha et al.|2004), ghost inflation ( Arkani-Hamed et 
al.|2 004) or trapped inflation ( [Green et al 



2009). Flattened- 
triangle configurations arise from initial conditions which 



deviate from the Bunch-Davies vacuum (Holman fc ToL 
|ley||2008l ). The limits for /jvl of the local shape from the 
WMAP5 analysis are 
(i.e. 2a) ( [Komatsu et al.||2009a l 



9 < fNT' < 111 at 95% confidence 



After horizon crossing, non-linearities in both the gravi- 
tational and hydrodynamical evolution of the baryon-photon 
fiuid prior to recombination ( Pyne fc Carroll|1996 Bartolo, 



Verde 2009 1 or the Sunyacv-Zcl'dovich effect ( [Goldberg fc 
Spergel[|1999 \, as well as foregrounds such as dust, galactic 



synchrotron radiation and unresolved point sources. Finally 
there are also instrumental effects, see e.g. [Donzelli et al.[ 
(2009). These effects contribute spurious non-Gaussian sig- 
nals, thus biasing the measurement of the primordial signal. 
It is therefore important to develop tools to isolate the pri- 
mordial signal from the contaminants. 

The separation of the CMB and foreground compo- 
nents exploits their differing spectral and spatial distribu- 
tion (Goldberg & Spergel 20081. It has been shown that 



foreground sources distributed in an anisotropic manner can 
mimic a non-Gaussian signal which could be taken for pri- 
mordial ( jCarvalho^ ,2009) . It has been suggested that the 
claimed detection of /jvl ( [Yadav fc Wandelt|2008[ ) might be 
due to residual foregrounds or instrumental contamination 



( Pietrobon et al.|2008{[Cabella et al.|2009[ [Smith, Senatore 
[fc Zal darriaga"2009') . Consequently it is important to de- 
velop estimators based on the bispectrum also able to mask 
out non-primordial parasite non-Gaussianities. 

In this paper we shall focus primarily on bispectral non- 
Gaussianity of the "local" type just defined, although most 
of the discussion generalizes straightforwardly to other tem- 



plates. See also Fergusson fc Shellard (2007 20091, whose 



authors are working on another bispectral estimator pro- 
gram, specifically aimed at dealing with theoretical primor- 
dial bispectra of general momentum dependence. In fact, 
their expansion of a general primordial bispectrum in terms 
of separable eigenfunctions is one way to easily extend our 
work beyond the "local" type of primordial non-Gaussianity. 
Unlike much of the literature, where the emphasis is on de- 
veloping an "optimal" estimator, where all the data are re- 



duced to a single number ( 


Heavens 


1998 


Komatsu et al. 


2005 


[Yadav et al.[[2008| [Creminelli, Senatore fc Zaldariaga 


2007 


1, the emphasis in our paper is placed on developing a 



more fiexible approach where the data can be divided in 
many ways and under which the effects of parasite non- 
Gaussianities can be masked in a well-defined way. See also 
Curto, Martinez-Gonzalez fc Barreiro' (2009);'Lan fc Marin- 
ucci[([2 008|);| Rudjord et aTTp jOO); Hika gc (2008);|Munshi fc 
Heavens (2009); [Munshi et al. (2009); Santos et al.[ ([2003[); 



Smith, Senatore fc Zaldarriaga, ( ,2009) for other estimators, 
many not based on the bispectrum but on wavelets, needlets, 
or Minkowski functionals. 

The paper is organized as follows. In section 2 we es- 
tablish some notation and discuss some general properties 
of estimators. A simplified and intuitive discussion of the 
nature and distribution of statistical weight of the "local" 
bispectral signal is given and it is shown how to modify 
templates in order to best mask parasite contaminants. Sec- 
tion 3 presents a quantitative discussion of the shape of the 
"local" signal and in particular the acoustic oscillations are 
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examined quantitatively. In section 4 the signal-to-noise for 
acoustic oscillations in the vertex region (where £i ^ ^2, ^2) 
is examined. Section 5 gives the computational details of 
a binned bispectral estimator that can be used for a gen- 
eral template and presents tests validating this estimator. 
Finally in section 6 we present some concluding remarks. 



2 BASIC FORMALISM 
2.1 Reduced bispectrum 

The temperature fluctuations in the CMB T{Q) = 
ATiyi) /To can be decomposed into spherical harmonics ac- 
cording to 

T(n) = ^a^™y,„(n) (4) 

with the harmonic coefficients aim given by 

aim = / dnT{Q.)Yim{^). (5) 



The bispectrum on the celestial sphere consists of cubic com- 
binations of the spherical harmonic coefficients of the form 



(6) 



whose expectation values may be calculated for a given the- 
ory. However under the assumption of statistical isotropy, 
these expectation values are not independent and can be re- 
duced to quantities depending only on £2, and £3. In the 
sequel, we shall not always explicitly indicate expectation 
values, because in many cases the equations can be inter- 
preted either with expectation values or as relating statistics 
of a particular realization. 

We may define a manifestly rotationally-invariant re- 
duced bispectrum in terms of integrals of triple products of 
maximally filtered maps so that 



h,e,e, = / dCl Tt,{Cl) Tt,{Cl) Te,{Cl), 



(7) 

where the maximally filtered map is defined as 

Ti{p.) = atmYim{(l). (8) 

m=-£ 

Using the expression for the Gaunt integral 

STXtT" = I d^Yl,m,(n)Y,,mmYl,m,m (9) 



(2^i-H)(2f2 + l)(2£3 + l) 
47r 

£2 £3\ f £1 £2 £3 



we obtain 



J \rni 1712 rn3 J ' 



(2£i 4-1) (2^2 + 1) (2^3 + 1) f£i £2 £3 
47r VO 

h £2 £3 



El -1 1 j^r, 

\mi m2 ms/ ^i'2<:3 

mi ,m2 ,m3 



^2 <-3 \ jmxm^m-i j-^g-j 

As a consequence of the Wigner-Eckart theorem, fe^^^g™^ is 



proportional to 



( £1 



£2 



\mi 7712 "l3 

with the 3j-symbol identity 



Using this fact combined 



E 



mi mo TTiQ 



£1 £2 £3\f £1 £2 £3 
nil m2 m3 J \mi m2 m3 



(11) 



which holds whenever £i, £2, £3 satisfy the triangle inequality 
and the parity condition (—1)^1+^2+^3 _ (-(-1), we find that 



1 mi m2 

°«i<2<:3 



Ivr 



(2£i + l) (2^2 + 1) (2^3 + 1) 



(12) 



£1 £2 £3 




1 

mi 



£2 
m2 



f.3 

m3 



Appealing to the flat sky approximation, we may inter- 
pret I?! as the expectation value for a single triangle with 
sides of dimension £1, £2, and £3 multiplied by the number of 
possible such triangles on the celestial sphere N/s{£i, £2, £3)- 
The number of triangles A'^a {£1 ,£2, £3) may be calculated by 
considering the variance of be-^^i^^t^ in a Gaussian theory, 
which should be given by 



Yarlbi^e^ta] = A^a (^1 , £2, £3) ci^^Ci^ci.:^ g, 



(13) 



which may be regarded as the formal definition of 
N/s{£i,£2,£3)- Here g is a combinatorial factor, equal to 1, 
2, or 6 depending on whether one, two, or three of the £ are 
equal, respectively. Using the above formulae, we obtain 



iVA(£l,£2,£3) 



(2£i + l)(2£2 + l)(24 + l) M £2 £3 



Itt 







(14) 

We now obtain a second useful reduced bispectral quan- 
tity, which is often called the reduced bispectrum in the lit- 
erature. 



NM£i,e2,£3) 



(15) 



4tt 



{2£i + 1)(2£2 + 1)(24 + 1) 



£2 




E 



£1 £2 £3 

mi 7712 ni3 



T mim2m3 
°fl<'2^3 ■ 



This definition is the direct analogue of b((.i,£2,£.3) in the 
flat-sky approximation. 



2.2 Natural inner product and ideal estimator 

Under the assumptions of full sky coverage, isotropic instru- 
ment noise characterized by a power spectrum rii (where 
the beam proflle has been accounted for by augmenting the 
noise) and the absence of foregrounds, the following opti- 
mal estimator may be constructed for /jvl (assumed to be 



^ For example, for uncorrelated (i.e. white) detector noise and a 
Gaussian beam profile 

m = no exp [£{£ + 1)6 f 

whm 

'/(81n2)] . 
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small) 



/"ideal 
NL 



(16) 



Here the natural bispectral inner product (relative to a par 
ticular experiment characterized by ne) is defined as 



(17) 



where the variance Var[6^J'|2£3] is calculated according to 
eqn. (13 1 except that ci has been replaced by {ct + ni). [In 
the remainder of this paper Var[b] will always include the 
noise ni.\ In eqn. (16 1 6°'"' represents the observed bispec- 
trum, which here is rescaled to correct for the attenuation 
from the finite width beam profile. Note that we are here 
using the Gaussian assumption to calculate the variance, 
which is only correct in the limit of weak non-Gaussianity. 
In fact , as we will see in section|5.3|and as has been pointed 
out in Creminelli, Senatore & Zaldariaga (20071; Liguori et 



al. (20071, if the value of Jnl corresponds to a detection at 



several sigma for a given experiment, non- Gaussian correc- 
tions to the variance become important. 

The estimator above can be rewritten in the form 



/"ideal 
NL 



5Z ™^1*2«3 
Il<ii2^t3 



lobs 
O1I1I2I3 



(18) 



which makes apparent the fact that the optimal estima- 
tor is an inverse variance weighted linear combination of 
the estimators that would be obtained for each independent 
(^1,^2,^3) triplet. To leading order, estimators correspond- 
ing to distinct triplets are uncorrelated. The inner product 



(6*", 6'") defined in (171 is equal to the total x or {S/Ny 



for detecting a bispectral template h^^i^i^ for an experiment 
whose noise and beam profile are defined by . 

The estimator presented here assumes uniform sky cov- 
erage. However many scanning strategies, and in particu- 
lar the one adopted for Planck, cover some parts of the 
sky more densely than others. For the present estimator in 
the form described here, such non-uniform noise will lead 
to an increase in variance but without any bias. Because of 
the non-uniformity of the essentially white instrument noise, 
products of map pairs, particularly those of large wavenum- 
ber, will be correlated with the scan pattern, which includes 
power at low wavenumber. This means that correlations 
with a third map of low wavenumber (i.e., the triplet com- 
binations contributing predominantly to the total available 
signal-to-noise squared for local bispectral non-Gaussianity) 
will include a product of the low-^ projection of the scan 
pattern with the low-^ Gaussian component of the primor- 
dial anisotropy, the expectation value of which is zero be- 
cause the scan pattern is not correlated with the primor- 
dial anisotropy. This effect can be removed by calculating 
the expectation value of the product arising from the scan 
anisotropy and subtracting it. The effect of non-uniform sky 
coverage is not necessarily all harmful, because deeper cover- 
age of some parts of the sky provides additional information 
which in principle could be exploited to obtain less noisy 
correlations between the low-^ primordial anisotropy, on the 
one hand, and the local high-£ power, on the other hand. 



The above problem and essentially the same strategy for 
mitigating it was first noted and implemented in |Creminelli| 
([2OO6I). 



et al. 



2.3 Qualitative nature of signal 

To understand intuitively where the statistical information 
for the non- Gaussian signal predicted from inflationary mod- 
els is situated, it is instructive to express (S/N)^ as a sim- 
ple integral derived under the following simplifying assump- 
tions. We employ the flat sky approximation, ignore the dis- 
creteness of the multipole numbers, and set all combinatorial 
factors to one. We further assume that eg ~ £~'^ (in other 
words, a power spectrum devoid of acoustic oscillations, Silk 
damping, and a flnite width for the last scattering surface). 
We assume bi^i^i^ ~ (^r^^2 ^ + ^2 ^^3 ^ + ^3 ^^7^) > which 
follows from exactly analogous assumptions. Taking into ac- 
count only cosmic variance, we obtain the triple integral 



^2 



(19) 



The integral obviously should be cut off both at small I (be- 
cause the quadrupole is the lowest accessible multipole) and 
at large I as well, because the finite resolution of the ex- 
periment becomes relevant and serves as a cut-off to render 
the integral finite. By simple power counting, one would ob- 



tain {S/Ny 



but considering elongated triangle 



configurations, with l-i <C ^2,^3, one obtains a logarithmic 
divergence factor, refining our estimate to become (see also 
Babich fc Zaldarriaga| (|2004[)) 



■ ^sky^max 



In 



(20) 



We suppressed the dependence on /ml which would add a 
factor fNL^- The above estimate, despite the crudeness of 
the approximations employed, highlights a number of qual- 
itative features of the non-Gaussian signal predicted from 
inflation. The presence of the logarithm emphasizes the im- 
portance of the coupling between the largest and smallest ac- 
cessible scales of the survey. The n^fcy^^^^ is proportional to 
the total number of useable modes (or equivalently pixels), 
which would suggest that most of the statistically signifi- 
cant information lies on scales at the limit of the resolution 
of the survey. If there were no logarithmic factor (in other 
words, no infrared divergence in the above continuum inte- 
gral), little information would be lost by dividing the full- 
sky map into a large number of smaller maps and forming 
a consolidated estimator by averaging the /jvl's estimated 
from each of these submaps. However doing so would entail 
throwing away the logarithmic enhancement factor, whose 
presence indicates the importance of the correlations be- 
tween the largest and smallest scales of the survey. If it were 
not for acoustic oscillations, which reverse the sign of the 
predicted bispectrum with respect to the simplistic model 
adopted above, one could construct a crude caricature of 
the ideal estimator by dividing the survey in two maps, so 
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that 

and considering the correlation 



(21) 



(22) 

whose expectation value would vanish in the Gaussian the- 
ory, and which measures to what extent the high-^ power 
is modulated by the low-E anisotropy pattern. We shall see 
later on that this estimator is far closer to the actual esti- 
mator than one might believe at first sight. 

2.4 General linear /jvl estimators 



The inner product defined in (17 1 is useful for calculating 



how much information is lost by employing a non-ideal tem- 
plate to estimate /jvl rather than the exact theoretical tem- 



plate appearing in (161. Suppose that we use the template 
j^tempi form the estimator (where the factor in the de- 
nominator provides the normalization needed to render the 
estimator unbiased) 



non — ideal 



tempi Lobs 



^ ' JJVI 



(23) 



In principle any template with a non- vanishing overlap with 
the theoretical prediction could be used. The lack of opti- 
mality is characterized by the ratio of variances 

,th 



TT- / rnon — ideal\ 

Var(/«£'") 



/■Lth ■Lth \ 



(24) 



indicating by what factor the variance is increased relative to 
the minimum variance estimator. Here 9 is the angle between 
the directions determined by the two templates in template 
space. 

There are several reasons why using a non-ideal estima- 
tor could be attractive. Firstly, to compute the ideal esti- 
mator exactly would require an inordinate amount of com- 
puter time, most of which would be of marginal value as 
the A£ = 1 spectral resolution is not needed. Secondly, in 
the presence of contaminants the above weighting scheme 
is no longer ideal anyway and another weighting must be 
employed. Thirdly, let us assume that a detection is made 
by means of a single number. Then we would want to divide 
data in many different ways to make sure that the detected 
signal is not of another type with fortuitous overlap with 
our favorite theory. For a detection with significant signal- 
to-noise showing that subdividing ^i£2^3 space in different 
ways always gives a consistent Jn l would lend credence to a 
claimed detection. These points will be investigated in more 
detail below. 

The ideal estimator is time consuming to calculate. If 
we cut off the sum in harmonic space at, say, H ~ 10^, we 
need to consider ~ 10^ reduced bispectrum terms (actu- 
ally the true number is about a factor 10 smaller because 
of the triangle inequality) and the calculation of a typical 
observed bispectrum coefficient would involve ~ 10® terms. 



KSW(Komatsu et al. 20051 discovered a clever algorithm 
for calculating the above estimator, but this method is re- 
stricted to producing just the single number /jvl, and has 



fixed weights. To go beyond these restrictions another way 
to speed up the computation is desirable. 

Obviously the data analysis challenge just presented is 
exaggerated, because the limit on the spectral resolution is 
without relation to physics underlying the CMB, and the 
scales over which interesting variations in the predicted bis- 
pectrum may be expected to occur are characterized by the 
period of the acoustic oscillations « 200, and the damp- 
ing scale Idamp ~ 400, which is due both to Silk damping 
and the finite width of the last scattering surface. It is not 
necessary to compute the bispectrum for every ^-triplet. One 
can instead bin in harmonic space to drastically reduce the 
computational effort. 



2.5 Binning 

We define a coarse-grained, or binned, estimator by parti- 
tioning into a sequence of M bins 



A+1 



11 



(25) 



where A, B = (1, . . . , M) and £i < £2 < . . . < Im+i (the last 
bin includes Im+i)- We define the binned bispectrum as 

hABC = X] 5Z 5Z b'^A^Btc (26) 

and an analogous expression for the binned variance 
Var[6^Bc]- 

The increase of the variance due to binning can be quan- 
tified using the cos^ 9 defined above in ( 24 1 . If binning is 



used, it corresponds to the following bispectral template: 



Var[&fi 



Var|6^ecJ 



(27) 



where (I^,Xb,Ic) in this expression is that particular bin- 
triplet that contains the mode-triplet {t.i,t.2,(-i)- Results for 
different binnings are presented later in the paper. Note that 



{b ,b ) = {b 



'.nned 7 binned\ 

,b ). 



(28) 



2.6 Dealing with bispectral contaminants 



In the presence of contaminants the above weighting scheme 
is no longer ideal and another weighting must be employed 
(see also ( Smith fc Zaldarriaga|[2006 l). This situation may 
be illustrated by the following simplified model, where the 
observed bispectrum may be decomposed into the sum 



7 ohs J- j th 

b = /jvlo 



(29) 



Here all the b's have three subscripts £1, £2, £3, which have 
been suppressed to render the notation more compact, and 
fe'^""' is a contaminant bispectral component (for example 
due to incompletely subtracted foregrounds or instrumental 
effects) and jg the noise from a particular realization 

of the Gaussian theory. In order to modify our estimator to 
mask the effect of the contaminant, we must characterize 
the statistical properties of the contaminant bispectrum. It 
is reasonable to assume that 



(6"°"*) = 0, 



(30) 
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because any known bias can be removed by direct subtrac- 
tion from the observed bispectral signal. The simplest statis- 
tical description for fo'^""* would be a Gaussian, which would 
be completely characterized by the correlation function 



^cont ncont rcont \ 



(31) 



We compute the likelihood ignoring constant factors not de- 
pending on Jml and integrating out the unknown contami- 
nant 6"°"* 



C 



{Vb' 
X exp 



)exp 



2 



1/7 ohs icont r j i 

-^{b -b - fNhb 



thsT 



- Jnlo ] 



(32) 



which may be evaluated to (up to an irrelevant factor) 
1 



exp 



obs r I th\T / r~^G , /^cont \ — 1 / 7 ohs r rth\ 

Jnlo ) (C + C ) (o — jNLb ) 



which we would like to put into the form 



' exp 



1(/a 



Jnl) 



/jVZ, 



(33) 



(34) 



We obtain the following optimal estimator for this modified 
situation 



fml _ b 

Jnl 



and the "information" is given by 



-3!— = b"'^(C'^ + c'=°"t)-i6*h^ 



(35) 



(36) 



3 DESCRIPTION OF THE THEORETICAL 
BISPECTRUM 

We now turn to examining more quantitatively the char- 
acter of the predicted reduced bispectrum b(^i,^2,^3) that 
results from the "local" bispectral anisotropy as a result of 
the ansatz ([sjl. Plots of the primordial bispectrum for var 



ious shapes can be found in Fergusson & Shellard (2007 



2009), both 2D slices and 3D, but here we focus in detail 



on the local shape. All simulations in the rest of this paper 
assume cosmological parameters from the WMAP five-year 
best fit |Komatsu et al.| (|2009a[) using WMAP data only. In 



Sect. |2.3| we examined the signal-to-noise ratio making the 
most brutal approximation so that all relevant quantities 
could be expressed in terms of simple analytic expressions. 
In particular, under the assumptions made there, we found 
that 



(37) 



We would like to study the shape of the local bispectrum 
more carefully in order to understand the deviations from 
this simple functional form, which moreover for plotting is 
extremely useful because it converts the bispectrum to a 




Figure 1. We plot the reduced CMB TTT bispectrum rendered 
dimensionless in a scale- free way according to the function B^-^i^e^ 
on 12 sections of constant {£i + 12 + ^3)/3 corresponding to 200, 
300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300. The 
color scale (explained in the text) ranges from —14.0 to +14.0 
and values outside of this range are clipped. 



quantity that is comparatively rather constant. Therefore it 
is useful to define and plot the dimensionless quantity 



B 



(38) 

where for purposes of normalization /nl ~ 1 has been as- 
sumed. The normalization constant a has been chosen such 
that cj^ — a£{£ -I- 1) in the Sachs- Wolfe approximation 
AT/T = —(1/3)$ where $ is the primordial Newtonian 
potential. 

Another choice (leading to a quantity that is even more 
nearly constant) is 



^^i^2^3 = , " ' , . (39) 

Here using the actual power spectrum in the denomina- 
tor factors out the overall decay (i.e., due to Silk damp- 
ing and the finite width of the visibility surface). However, 
this quantity has the disadvantage that the oscillations in 
B become intertwined with the oscillations of the two-point 
power spectrum. 

The panels of Fig. [ijplot Be-^e^ts , showing a sequence of 
sections of constant {ii + £2 + £3)- Fig. |2] shows the analo- 
gous plots for B. The xy-coordinates in the plot correspond 
to the orthogonal transverse coordinates £±a ~ {£1 ~£2)/V^ 
and £±1, — (£1 + £2 — 2£s)/\/6. The function value is indi- 
cated by height and color with the color scale ranging from 
-14.0 to -1-14.0 for B and from -3.0 to -1-3.0 for B. Black -s> 
blue — ^ green indicates negative values and red — > orange — >■ 
yellow indicates positive values, with the zero point situated 
between green and red. 

The triangle inequality (|^i— ^2! ^ £3 ^ £i+£2 and cyclic 
permutations) reduces to the following three conditions 

Li=£2+£3-£i ^ 0, 
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Figure 2. Same as in Fig. ^except that we instead plot 3*1^2*3' 
normalized with the actual CMB temperature power spectrum 
rather than the scale free one. Here the color scale (same as above) 
ranges from —3.0 to +3.0. 



L2 = 4 + ^1 - ^2 > 0, 



(40) 



Because the sum (£\ + £2 + £3) has to be even (parity con- 
dition), the Li only take even values. The inverse relations 
are 

L2 + is „ ii + L3 ii + L2 



£1 



(41) 



2 ' 2 ' 2 

If all non-negative values of £1, £2, £3 were allowed, the 
region allowed by the triangle inequality would comprise 
the infinite triangular pyramid in the principal octant with 
£i,£2,£3 > whose edges are the rays = ^2, ^3 = 
and the two other cyclic permutations thereof, and the faces 
would correspond to Li = (constant). The extremely elon- 
gated triangles with one £ much smaller than the other two 
are situated in the immediate neighborhood of these edges. 
It is these configurations that give rise to the logarithmic 
divergence previously mentioned and contribute the bulk of 
the statistical weight, in a way which we will characterize 
more quantitatively later on. Because the lowest observable 
multipole starts at ^ = 2 the edges are slightly truncated by 
the intersecting planes {£1 ^ 2, . . .). Fig. [ijshows 16 sections 
of constant {£i+£2+£3)- The most negative values of B occur 
near the edges. In the central region we observe a sequence 
of acoustic oscillations. First, we observe a single maximum 
situation at the center of the triangle (corresponding to an 
exactly equilateral configuration). The center oscillates back 
down but three new maxima appear situated between the 
centre and the vertices. The process repeats itself, so that 
the number of local maxima progresses according to the se- 
quence 1, 3 = 1-^2, 6 = H-2-1-3, 10 = 1-^2 + 3-^4, ... [An 
animation of a sequence more finely sampled is available at 
http://www.apc.univ-paris7.fr/~bucher/bispectrum,] 

In Figs. [3] and |4] we show isosurfaces of the function 
Beii2i3 at values -10, -4, -2, +2, -1-4, -1-10. The domain is 
an infinite triangular pyramid in the first octant and in the 
first panel all the isosurfaces for this most negative value 



terminate on the boundary. The long curved lens-like struc- 
tures are situated along the edges and track the acoustic 
oscillations of the two-point power spectrum, as discussed 
in more detail below. As the threshold isovalue is increased, 
new bubbles appear and coalesce. The null isosurface is not 
shown because it has a regular Swiss cheese appearance, so 
that it is impossible to see inside. 

The TTT reduced bispectrum is given by |Komatsu fc| 



Spergell ( |2001 1 



ki^dki 



k2 dk2 
Jo 



X Ae^ki) Ai2{k2) Ai.,{k3) 



P{kl)P{k2) + P{k2)P{k3) + P(fc3)P(fcl) 



X/ r dr ji^{kir) jt2{k2r) jtsiksr), (42) 
/o 

where the primordial power spectrum for the gravitational 
potential <1> is defined by 

{$(k)-l'*(k')) = i2nf5^{k - k')P(fc) (43) 

and the local ansatz ([3| is assumed. Ai{k) is the CMB trans- 
fer function (computed with a Boltzmann solver, here using 
CAMB). At the same time the expectation value for the 
CMB power spectrum is given by 

2 



ce = 



k^dkP{k)Ai{k). 



(44) 



If we assume a simple power spectrum P{k) (e.g., exact 
scale invariance or a homogeneous power law), its contri- 
bution to the reduced bispectrum is featureless and does 
not single out any special scales. However the CMB transfer 
function Ai{k) contains interesting physics characterized by 
a variety of scales (e.g., the scale of matter-radiation equal- 
ity, horizon size at last scattering, the shape of the visibility 
function, etc.) and it is these features, which are presumably 
well understood, that should be helpful in distinguishing a 
primordial bispectrum from other contaminants. 

We re-examine some of the arguments previously put 
forth by exploiting the approximate integral for the statis- 
tical information in eqn. ( |19[ ) in light of the exact local bis- 
pectrum and exact inner product for a particular experiment 
(here taken to be the PLANCK instrument using an uncon- 
taminated CMB map obtained by combining in quadrature 
the 100, 143, and 217 GHz channels with the specifications 



published in the PLANCK bluebook (Planck Collaboration 
[2006] )). 

In order to visualize how the statistical weight of the 
signal is distributed in £\£2£3 space, it is useful to make 



various cuts and integrate x fj^ 



taking into account these 



cuts. (We take /jvl = 50 as that is roughly the order of 
magnitude of the central value of the WMAP5 results.) The 
simplest cuts one can make is to restrict the multipole values 



€ £ € 



The resulting x 



to lie in a range £rmn 
function of £m.in and £max is shown in Fig. [S] We observe 
that almost all (« 2/3) the statistical weight comes from 
the coupling of the lowest £ {£ <y 10) with the largest £ 
(^ 500). 

It is useful to adopt the following parameterization un- 
der which there is one dimensionful quantity and a number 
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Figure 5. Wc indicate the total x for ^ detection of Jml = 50 
for the Planck experiment taking into account the three frequency 
channels 100, 143, and 217 GHz. In the upper plot is shown as 
a function of imin and imax when only multipoles in the range 
imin ^ ^ ^ imax are taken into account. In the lower panel, a 
one-dimensional cut of the dependence on imin (while £ 

max is 

held at its maximum value of 3000) is shown. 



of redundant dimensionless quantities (of which only two 
are linearly independent). For our dimensionful quantity we 
choose isum = ^1+^2+^3 ~ L1+L2 + L3. The dimensionless 
shape of the triangle in harmonic space is characterized by 



Xi 



Li/is 



X2 



X3 = 



(45) 



where Li , L2 , and La are as defined in eqn. ( 40 1 . The quadru- 
plet {laum,Xi,X2,Xz) provides a redundant set of coordi- 
nates (satisfying xi + X2 + xz = 1) and the triangle inequal- 
ity is satisfied if and only if ^ a;i ^ 1, 5S a;2 5S 1, and 
^ 2:3 5C 1. The three edges of the allowed region are defined 
by (L, 1,0,0) and (L, 0,1,0) and (L, 0,0,1) where L ^ 6. 

Fig.jejshows the weights {h^'^ f /V&r[b°*"] (the sum over 
£1,^2, is of which corresponds to the x^) for three different 
slices with iaum fixed at 1000, 2000, 3000, respectively. Since 
the absolute value of the weights on such a slice is not very 
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Figure 6. To illustrate how the statistical weight for a detec- 
tion of /jvL is distributed, we plot the weights (6*'')^/Var[f)°''^] 
(for the same three Planck channels as in the previous figure) 
on three representative slices: £sum — 1000 (top), isum — 2000 
(middle), and isum = 3000 (bottom). The normalization is such 
that unity corresponds to the average value on the slice. In the 
central regions the grey contours correspond to sign changes of 
the bispectrum. Note that the condition ^1 ^ ^2 ^ is means that 
only 1/6 of the full triangle is relevant, the rest consisting of 5 
identical reflected copies. 
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meaningful, they have been divided by the average value on 
the slice. In other words, values smaller than unity corre- 
spond to regions with lower than average weight, and values 
larger than unity to regions with higher than average weight. 
The preponderance of the extreme configurations around the 
vertices is manifest in this figure, although their integrated 
contribution is not so easy to see. 

We now define two families of cuts in order to quantify 
how much the triangles of the central region contribute to 
the total integrated /nl signal. Firstly we may define a fam- 
ily of expanding concentric equilateral triangles parameter- 
ized by A according to the relation A ^ a;i, a;2, 2:3 ^ (1 — 2A) 
with ^ A ^ 1/3. This is one way to cut out nearly flattened 
or coUinear triangles. Secondly, we may cut away at the cor- 
ners by requiring that i'l/^sum, ^2/4tim, > 7 where 
I (-sura ^ 7 5S 1/3 SO that the triangular region becomes 
hexagonal. The results of the integrated for /jvl = 50 for 
these two families of cuts is indicated in Fig. [7] In the second 
panel the integrated contributions for positive and negative 
values of the bispectrum are examined separately. We ob- 
serve that the exploitable signal is heavily concentrated in 
the corners and very little statistically useful signal is present 
where the bispectrum is actually changing sign. 



4 PROSPECTS FOR DETECTING ACOUSTIC 
OSCILLATIONS FROM VERTEX 
CONFIGURATIONS 

In the previous section we investigated the imprint of physics 
occurring around the time of recombination on the bispec- 
trum, which intuitively may be characterized by departures 
from a constant function of the bispectral ratio Bi^t^i^ de- 



fined in ( |38[ ). We observed a rich pattern of acoustic oscilla- 
tions in the a;ia;2a;3-plane as I varies. However unfortunately, 
the statistical weight of this central region is negligible, at 
least for an experiment like PLANCK given the limits on 
/jvz, set by WMAP. For PLANCK there is simply not enough 
signal-to-noise in this region, as we have seen, although there 
may be some hope with a next-generation satellite experi- 
ment such as some of the most ambitious proposals for a 
B-polarization satellite. 

The remaining window for detecting bispectral oscil- 
lations resides in the ^-dependence of b^j^jfj in extreme 
vertex configurations. These configurations are defined by 
^1 ^ ^2,i'3, which automatically means that (2 ~ be- 
cause of the triangle inequality. Hence we define Ismail = (.1 

and harge = (^2 + ^3)/2 

this limit b. 



2 ~ ^3- Using definition (391, in 
can be written as 



2B{lsmall) ce^ 



(46) 



Here we have used the fact (shown numerically) that for 
these squeezed triangle configurations B is practically con- 
stant as a function of liarge- We conclude that in the 
squeezed limit the bispectrum is proportional to the power 
spectrum as a function of iiarge- This is shown explicitly 
in Fig. [s] [Alternatively, in a description that keeps the 3 
different multipoles explicitly, one can use the Li defined 
in (|40|. Then the squeezed configuration corresponds with 






Figure 7. The upper panel shows the effect on the total when 
one cuts away extreme triangles where one side is much shorter 
than the others. In the middle panel the contributions from pos- 
itive (blue, the curve that stays low) and negative (green) values 
of the bispectrum are shown separately. The bottom panel shows 
the effect of cutting away on the sides rather than the corners. The 
minimum values of 7 and A (defined in the text) correspond to 
taking all triangles into account, while for their maximum value 
(= 1/3) only equilateral triangles are taken into account. The 
same three Planck channels'^s 
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L2,Ls <C Li, the prefactor in (461 depends on L2 and L3 



only, and we have a power spectrum behaviour as a function 
ofLi/2.] 

To explain why B does not depend on iiarge in the 
squeezed limit, we start from the integral expression (421 
for the bispectrum b^^^j^g. In this limit the dominant con- 
tribution to the integral is proportional to 



dki 

'^i Jo Jo \k2 

X A<i(fci) Af,(A:2) Af3(fc3) 



X / r dr ji^{kir) ji2{k2r) ji^ikzr) 
lo 



fc3^ 



(47) 



where for clarity we insert an exactly scale invariant power 
spectrum, with P{k) ~ k~^, but all the arguments carry 
over to approximately scale invariant spectra. The dominant 
contribution to the radial integral in the last line occurs at 
r ~ a, where a is the radius of the surface of last scatter, and 
the dominant values of the radial wave number are situated 
around ki « £i/a. The first spherical Bessel function varies 
slowly while the the last two oscillate rapidly. This causes 
the integral to cancel except near the resonance, where k2 
is very nearly equal to fcs modulo a negligible offset of ap- 



proximately fci . We therefore approximate eqn. ( 47 1 by 



00 j7. 



k2^ k^-^ 



1^/ k2' 
'^i Jo Jo 

X Afi(fci) Ae^{k2) Ae^iks) 
S(k2 - k-i) r°° dr . 



k2 hi 



„ f°° dki ^ 
2/ ^A,,(fci 

Jo ''1 Jo 



(48) 



So we recover the fact that the bispectrum is proportional to 
the power spectrum as a function of iiarge, times a prefactor 
that only depends on Ismail- We have confirmed by numeri- 
cal integration of the exact expression that the factorization 
on the last line is a good approximation when ismaii is much 
less than the scales of the power spectrum acoustic oscilla- 
tions and the damping. This factorization property greatly 
simplifies constructing estimators in the vertex region by 
reducing the feature detection problem to essentially a one- 
dimensional problem. 

We construct filters to detect the presence of several 
peaks and troughs in the bispectrum at positions corre- 
sponding to the peaks and troughs of the two-point power 
spectrum, which is of course now well known experimen- 
tally. The mathematical problem of how best to construct 
these filters is not completely well-defined, so there is some 
arbitrariness in the choice of filters designed to pick up the 
various features. Nevertheless other filter choices can be ex- 
pected to yield results differing only slightly. We first de- 
scribe our filter templates as if they were acting on the 
two-point power spectrum and then indicate how to carry 
over these templates to the bispectrum problem, which is 
essentially a question of how to deal with the two additional 
transverse dimensions. 

We define the absence of features as the constancy of 



Ce = l{i + l)ct and construct two types of filters. The first 
takes a discretized second derivative of the power spectrum 
and is used to detect the peaks and troughs. We require this 
filter to give zero for a linear rise or fall in Ce. The second 
filter evaluates a discretized first derivative in order to detect 
a rise or fall in Ci. 

The first filter type evaluates an integral over some win- 
dow function times the second derivative of Ci, where the 
support of the window function for a peak, for example, 
would extend approximately from the trough on the left to 
the trough on the right. This is the optimal width because 
making the filter narrower would decrease the signal-to-noise 
by relying on too few noisy data points and a broader fil- 
ter would bring in cancellations from neighboring regions of 
opposite sign, likewise diminishing the signal-to-noise. One 
candidate for the second-derivative filter would use the sec- 
ond derivative of a Gaussian profile, but this choice has the 
drawback of an infinite support, which would have to be cut 
off. Here we use instead a piecewise linear filter profile. For 
the first-derivative filter we use a step function constantly 
positive to the left of a certain value of I and constantly neg- 
ative to the right. The central value, the width to the left, 
and the width to the right are the three parameters in this 
filter. The relative heights of the left-hand and right-hand 
sides are fixed by the constraint that the integral of the fil- 
ter should be zero. The profiles used are shown in Fig. [8] in 
relation to the predicted signal. 

For the implementation in ^1^2^3-space, we include only 
(■small ^ 30. Going any higher would only marginally im- 
prove the signal and would tend to smooth the power spec- 
trum. Then we group together triplets with the same Iiarge- 
The function f {Ismail) from the factorization in eqn. (48 1 



is found numerically and with each such group the best es- 
timator for llarge{ilarge + l)t'(^ small 1 ^large small 

) (the 

expectation value of which is C^,^^ ^) is constructed. The es- 
timators for each triplet are combined according to inverse 
variance weighting. Explicitly, the resulting estimator Ct,^^g^ 
is given by (using the more precise notation Li/2 instead of 

(large) 



c. 



Li/2 



iVi 



2 fti^^faiVA (^1,^2,^3)^1 



2)2^ LxlLx + l 



2^^1^30 



Var[&?^^»^,3] 



(49) 

where Li = £2 + £3 — £1 as before and the sum is over the 
two directions transverse to Li, i.e. L2 and L3, under the 
condition that £1 ^ 30. Here A'' is the normalization of the 
inverse variance weights, given by 



N 



> 2 



= E 

2<£i <30 



{2Bi,t2i3N^s{,£x,£2,£3)ci,] 



^(^ + 1)) Var[&°^^?^,3] 



(50) 



This estimator and its signal-to-noise ratio squared for the 
case of an /jvl = 50 detection are shown in Fig. [S] It repre- 
sents one way of going beyond the single number /jvl to de- 
scribe the bispectrum (see also Munshi fc Heavens| ( |2009 1). 

To look for specific features, the estimator Cti^^^^ is 
summed with respect to the filters J-{£iarge) described 
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Feature 


Center 




1 ' ^ V ) ^ 1-^ + 


(S/N)'^ 


1st peak 


210 


100 


170 


0.43 


1st trough 


400 


130 


130 


0.69 


2nd peak 


530 


100 


120 


0.19 


2nd trough 


660 


100 


120 


0.27 


3rd peak 


810 


130 


160 


0.70 


3rd trough 


1010 


140 


120 


1.08 


4th peak 


1140 


110 


120 


0.34 


4th trough 


1300 


130 


120 


0.53 


Drop A 


910 


580 


1770 


11.3 


Drop B 


1550 


1200 


760 


21.9 



Table 1. Signal-to- noise squared for detecting certain features 
in the edge bispectrum £i <^ ^2i^3 for /jvl = 50. Here we have 
cut with £i ^ 30. The first eight filters are second derivative 
filters designed to detect peaks and troughs. Drops A and B are 
broader filters representing a discretized approximation of the 
first derivative centered between the third peak and third trough 
and between the fifth peak and fifth trough, respectively. These 
first derivative filters give a much larger (S/N)'^ because of their 
extreme breadth. See Fig. |8] for the filter profiles. 



above. The resulting statistic S is given by 



S - 



Ll/2- 



The variance of S is given by 

Var[S] = ^.F^(^)iV(i^), 



(51) 



(52) 



The signal-to-noise ratios of this estimator in the case of 
/nl ~ 50 are given in Table [l] Once more we have assumed 
the three Planck channels at 100, 143, and 217 GHz, and 
we have used imax = 3000. For an arbitrary /nl the table 
entries in the last column should be multiplied by a factor 
of (/jvl/50)^ 



5 BINNED BISPECTRAL ESTIMATOR 
IMPLEMENTATION 

In this section we describe how to estimate f^L from a CMB 
map having been provided the details of the CMB experi- 
ment (instrument noise, beam profile, . . . ) and a theoretical 
model, from which we require the primordial power spec- 
trum P{k) and transfer functions Ai{k). As explained in a 
previous section, we need to evaluate the binned approxi- 



mation to the ideal estimator ( 16 1 or to some variant ( 23 1 



We can write the binned estimator as 

fNL = Yl 



WABc{b% 



, 7obs 
I "ABC 



with the binned weights given by 



WABc{b%^=l) = 



E 



'■VET 



(53) 



(54) 



arf&i 



Since we have changed the sum over ^ ^2 ^ £3 into a sum 
over ^ ^ ^3 ^ C, we have to divide the weights by 2 (or 6) in 
the case where 2 (or 3) bins are equal, to avoid overcounting. 




— C I estimator 

— its cumulative (S/N)"2 

— - filter F 





Figure 8. The bispectrum for triangles near the vertices of shape 
space where li <^ £2 ~ £3 follows the shape of the two-point 
power spectrum as a function of Li/2 £iarge, as shown in these 
plots of the estimator Cii^^g^ defined in | |49[ l (green curves in the 
top and middle panels). Its (S/N)'^ for an /jvl = 50 detection is 
shown in the bottom panel, while the other solid curve (blue) in 
the middle panel shows the cumulative (S/N)'^ (and the values on 
the vertical axis of this plot correspond to this curve) . In the upper 

panel the dashed curves ^f*^ (j-^K)^^!^ ^l&I^^ffj^^6ffii'^ (Wf-9(jb 
derivative filters used to derect the peaks and the troughs. The 
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Computing this estimator involves four distinct steps: 

(i) Computing the theoretical bispectrum b^'^i^i^; 

(ii) Deriving the binned weights wabc from it; 

(iii) Computing the binned bispectrum b'^^c from the 
map; 

(iv) Combining the above to yield /ml (or another esti- 
mator). 

The modular code we developed follows these 4 steps and the 



results from each step are saved separately. Hence step (iii 



for example, is completely independent of steps (i) and 



so that when trying out another mask, say, only the code 
corresponding to steps (iii) and (iv) has to be rerun. While 
we focus in this section on /jvl, the modular setup of the 
code makes it very easy to compute another bispectral quan- 
tity instead. For example, by replacing in step (ii) the /nl- 



weights defined above by the binned version of the weights 
defined in (511, we can compute the statistic S to detect 



acoustic bispectrum features. 



5.1 Implementation notes 

After getting the transfer functions from CAMB, to save 
computing time in step (i) eqn. ( 42 1 for the reduced bis- 



pectrum is not evaluated for all ^-triplets, but only for a 
representative (large) number. The other necessary values 
are determined by 3D linear interpolation. Our code also 
allows for the use of a 3D cubic interpolation scheme, but 
for the large number of exact "grid points" used we found 
no significant increase in accuracy compared to the linear 
scheme, which is faster and already very accurate. 

To compute the weight for a bin in step |(ii)| one can 
either evaluate all terms inside the bin and sum them ex- 
plicitly or, much faster, for large bins use a simple 3D in- 
tegration scheme (here based on tetrahedronisation of the 
bin-cubes). Both methods are implemented in our code. In 
the first case all the values of ce, Bi-^^i^ir^, and A^'a (^i, ^2, ^3) 
are pre-computed and saved to speed up any consequent cal- 
culations for the same cosmological model (at the cost of ad- 
ditional memory usage). However, because of the slow vari- 
ation of ^^1^2*3 5 the second (approximate) method is highly 
accurate and due to its higher speed and lower memory us- 
age generally to be preferred. 

Computing the binned observed bispectrum in step |(iii)] 
involves multiplying the observed aim with the inverse beam 
and window function operators, applying the masks (if any), 
and determining the maximally filtered maps TeiQ.) accord- 
ing to (|8|. However, it is in this step that binning is a nec- 
essary ingredient to reduce computing time, so instead of 
computing the maximally filtered maps for each individual 
I, they are only computed per bin, producing Ta{^). Finally 
the binned observed bispectrum b^lBC is computed using the 
binned analogue of eqn. Q. Healpix ( Gorski et al.|[2005 l is 
used to carry out the spherical harmonic transforms. 



Using the files with the results from steps (ii) and (iii) 
it is simple and quick to evaluate in step |(iv)] the sum ( |53[ ) 
and find the result for fnh- Its variance is also computed. 
In the case of multiple maps (different realizations of the 
same CMB sky) the result is computed in a number of dif- 
ferent ways to check consistency: by first combining the dif- 



ferent bins to find a result per map, and then combining 
the different maps; or by first combining the maps to find 
an averaged result per bin, and then combining the bins. 
Moreover, using the set of maps we can also compute an 
observed variance and use that in the weights instead of the 
theoretical variance. In particular this allows us to see the 
effects of deviations from Gaussianity in the variance, which 
become important when Jml is several standard deviations 
away from zero for the experiment under consideration. 



Steps (ii) and (iv) are fast (about 10 minutes and 10 



seconds, respectively, on a laptop with an Intel Core 2 Duo 
P8400 processor (using only one core)), while step (i) is 
slower (about 1.5 hours), all for Imax ~ 3000. If necessary, 
step |(i)| can be made faster by sampling the (at present prob- 
ably unnecessarily densely sampled) theoretical bispectrum 



more sparsely. Step (iii) is the slowest, taking about 45 min- 
utes for one map with irnax = 2000 and Healpix parameter 
?^sido = 1024. However, we have parallelized this step, so 
that multiple maps can be run in parallel, and further code 
optimizations are in progress. We run this code on a blade 
server with 16 quad-core Intel Xeon 5345 processors on 8 
blades (64 cores in total with 4 GB of memory per core; 
one core of these performs roughly the same as a core of the 
P8400 mentioned above), so that a set of 100 maps takes 
about 1.5 hours. 



5.2 Binning scheme 

To run our code we have to specify a binning scheme. We 
choose a divisional scheme, although other schemes are im- 
plemented as well. This means that we first split the in- 
terval [£min,£max] into two at the intermediate value lint ~ 
l^min -\- aoilmax — £min) , wherc the parameter ao lies between 
and 1. Next each of these intervals is split into two again 
using parameters 60 and 61, each rescaled to lie between 
and 1. This process can be repeated as many times as nec- 
essary. 

We can quantify how much the variance increases due 
to binning, compared with an ideal estima tor without bin- 
ning, using the parameter cos^ defined in 1 24 1. The results 



are summarized in Table [2] To find the maximum value of 
cos^ 6 we have varied the 7 parameters ao, 60, 61, co, ci, C2, 
and C3 (only 3 parameters for the case of 4 bins, and for 
the cases of 16 to 64 bins the remaining divisional param- 
eters have been fixed at 0.5) on a fixed grid and assumed 
(■max ~ 2000 (as well as the noise and beam characteristics 
of the Planck 100 GHz channel). We see that using 64 bins 
instead of all the 1999 ^-values in the interval [2, 2000] gives 
an immense reduction of computing time at the cost of only 
a 0.7% increase in the variance. Even smaller numbers of 
bins still give a quite acceptable increase of variance. 

As an explicit example, the maximum value of cos^ d 
for 64 bins was reached for the following parameter values: 
ao = bo = Co = 0.2, fei = C3 = 0.4, and ci = C2 = 0.5. This 
corresponds to the following bin separator values for £: (2, 
3, 5, 7, 9, 11, 13, 15, 17, 25, 33, 41, 49, 57, 65, 73, 81, 101, 
121, 141, 161, 181, 201, 221, 241, 261, 281, 301, 321, 341, 
361, 381, 401, 441, 481, 521, 561, 601, 641, 681, 721, 761, 
801, 841, 881, 921, 961, 1001, 1040, 1088, 1136, 1184, 1232, 
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Number bins Overlap (cos'^ 9) 



fNL Estimator Performance (no mask) 
100 simulations; lmax=2000; Plancl< channels 100, 143, 217 GHz 

40 r 



4 
8 
16 
32 

64 



0.608 
0.846 
0.925 
0.974 
0.993 



Table 2. Quality of estimator as a function of number of bins for 

Imax = 2000. 



1280, 1328, 1376, 1424, 1496, 1568, 1640, 1712, 1784, 1856, 
1928, 2000). 



5.3 Validation of estimator 

To test our binned estimator, we applied it to 100 CMB- 
only maps with £max ~ 2000 and nsidc = 1024, kindly pro- 
vided by Michele Liguori based on the WMAP5-only best-fit 



parameters (Liguori et al. 2007j ). To these maps we added 
the effects of uniform white noise and a Gaussian beam ac- 
cording to the specifications of the 100, 143, and 217 GHz 
channels of the Planck experiment (combined in quadrature) 
( Planck Collaboration|2006[ |. We have used 64 bins, with di- 
visional parameters such that cos^ 9 is maximal. 

The results can be found in Fig. [9] Without a mask, we 
recover for /nl = a value of —0.3 ±5. 3 and for /nl = 50 of 
50.3 ± 7.6. These results agree with those obtained with the 
KSW estimator. We also tested with a simple band mask of 
[—6.5,-1-6.5] degrees around the galactic equator and found 
—0.05 ± 6.3 and 49.0 ± 8.5, respectively. Here the masked 
area was filled with the average of the rest of the map and 
any remaining monopole was subtracted. A more detailed 
study of the influence of masks and how to treat them will 
be the subject of a future publication. 

The smooth curve in both plots is the predicted Gaus- 
sian distribution based on the calculated theoretical variance 
when nonlinear corrections are ignored. Its standard devi- 
ation is 5.4 without a mask and 5.7 with the band mask. 
The modest increase with the mask is due to the decreased 
sky fraction. We see that for /ml ~ 50, which corresponds 
to a detection of several sigma, the actual standard devi- 
ation starts deviating from the Gaussian one as nonlinear 
corrections become important. 



6 DISCUSSION 

The search for primordial non-Gaussianity is presently in 
its early stages. Non-Gaussianity of the local type is moti- 
vated by multiple-field inflation and if detected would have 
a substantial impact on our understanding of the physics 
of the early universe. At present we have only an upper 
limit with /^£"' in the interval [-9, 111] at 95% confidence 



according to the official WMAP5 analysis (Komatsu et al 
|2009a[ ), which constitutes a hint of a detection. A detection 
has been claimed in another analysis of the WMAP data by 



Yadav & Wandelt ( 2008 1 , and several other alternative anal- 




fNL Estimator Performance (mask [-6. 5, +6. 5] degrees) 
100 simulations; lmax=2000; Planck channels 100, 143, 217 GHz 





Sigma = 


5.73 




-0.05 +/- 


6.28 




49.0 +/- 


8.54 




20 40 
fNL 



Figure 9. The performance of the binned estimator described in 
the text is demonstrated on two sets of mock maps generated by 
Monte Carlo, one with /jvi = and another with /jvl = 50. In 
the upper panel the entire sky is used whereas in the lower panel a 
mask of ±6.5° about the galactic equator has been used, resulting 
in a modest increase in estimator variance. Gaussian fits to the 
theoretical mean and variance are indicated in the histograms. 



yses, such as using needlets (Lan & Marinucci 2008 Rud- 
jord et aL]|2009 1 or spherical Mexican hat wavelets (Curto, 



Martinez-Gonzalez fc Barreiro||2009|, ha ve been published 
See also|Hik agc ( 2 008[) ; |Creminelli et a.l] (|2006[ [2007| ; |Sena- 
tore, Smith & Zaldarriaga'f 2009); 'Smith, Senatore & Zaldar- 
riaga| ( |2009[ ). These analyses differ in detail but are broadly 
consistent. Note that a confirmed detection of a primordial 
/nl of order 10 would rule out the simplest single-field slow- 
roll inflation models. 

The PLANCK space mission promises to bring about 
approximately an order of magnitude improvement in sensi- 
tivity to fNL, so if the present hint of a signal is confirmed, 
we will have an abundance of signal-to-noise at our disposal, 
allowing the data to be divided in many ways in order to con- 
firm the predicted shape for the bispectrum and to exclude 
other non-primordial sources of parasite bispectral power. 

In this paper, we demonstrated how a binned bispec- 
trum estimator can be useful to evaluate the perfect tem- 
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plate with a modest number of bins, so that applying the 
estimator is computationally extremely efficient. Compared 
to other implementations of the optimal estimator, which are 
likewise fast, the present estimator has the advantage that 
it can be used to apply other templates, either for /jvl or to 
determine other bispectral characteristics. It will be impor- 
tant to develop models of the expected pattern of bispec- 
tral non-Gaussianity arising from foregrounds, instrumental 
artifacts, and foreground subtraction residuals. Even if very 
approximate, such templates will be invaluable for excluding 
non-primordial explanations for a possible future detection. 
We showed how using a binned estimator one should mod- 
ify the template for primordial non-Gaussianity of the local 
type to optimally mask the contaminants. 

We also investigated how the acoustic oscillations man- 
ifest themselves in the CMB bispectrum from local primor- 
dial non-Gaussianity and in particular evaluated the signal- 
to-noise for detecting various features. Although a rich pat- 
tern of non-trivial acoustic oscillations occurs in the central 
region of triangle shape space, we found that these oscil- 
lations contribute negligibly to the total available signal- 
to-noise. Ifowever, the oscillations having the same shape as 
the power spectrum near the vertices of triangle shape space 
offer a much larger signal. We also investigated broader 
features due to the damping, and these may be seen with 
PLANCK if /jvl is large enough. An observation of such fea- 
tures in the bispectrum matching the shape of the two-point 
power spectrum would constitute a convincing confirmation 
of a first /jvl detection. The absence of a signal in other re- 
gions where none is expected also provides a powerful check 
against foreground and instrumental contamination. 
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